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ABSTRACT 

While tens or hundreds of stellar-remnant black holes are expected to form in globular 
star clusters, it is still unclear how many of those will be retained upon formation, and 
how many will be ejected through subsequent dynamical interactions. No such black 
holes have been found in any Milky Way globular cluster until the recent discovery of 
stellar-mass black holes in the globular cluster M22 (NGC 6656) with now an estimated 
population of 5 — 100 black holes. We present a direct TV-body model of a star cluster 
of the same absolute and dynamical age as M22. Imposing an initial retention fraction 
of ss 10% for black holes, 16 stellar-remnant black holes are retained at a cluster age 
of 12 Gyr, in agreement with the estimate for M22. Of those 16 BHs, two are in a 
binary system with a main sequence star each while also one pure black hole binary is 
present. We argue that multiple black holes can be present in any Milky Way cluster 
with an extended core radius, such as M22 or the model presented here. 

Key words: globular clusters: general - galaxies: star clusters: general - stars: mass- 
loss - methods: TV-body simulations 



1 INTRODUCTION 

Stellar-mass black holes (BHs) are formed at the endpoints 
of stellar evolution of very massive stars. Depending on the 
chemical composition or metallicity of a star, a supernova 
resulting in either a neutron star or a black hol e is the ul- 
timat e fate for most stars above ~ 7 — 8 Mq l|Pols et al.l 
Il998h . 

Tens or hundreds of such stellar-remnant BHs are ex- 
pected to form in globular clusters (GCs), however at the 
current stage it is still not entirely clear how many will re- 
ceive a velocity kick at formation and get ejected imme- 
diately, and how many BHs might be removed during the 
subsequent dyamical evolution of the cluster. While sev- 
eral extragalactic GCs containing BHs a re known at the 
current stage (e.g. iMaccarone et all 1201 lT ). none had been 
confirmed in a Milky Way GC until the rece nt discovery of 
two s uch stellar- mass BHs in M22 (NGC 6656, IStrader et all 
l2012h . Since only BHs currently undergoing observable ac- 
cretion can be dete cted via radio or X-ray emission and 
IStrader et ail (|2012T ) estimate that 2 - 40% of BHs are ex- 
pected to become members of binary systems with observ- 
able accretion over 10 Gyr, it is possible that a total popu- 
lation of w 5 - 100 BHs exists in M22. 

Theoretical predictions in the past have lead to the as- 
sumption that well before a cluster age of 12 Gyr all but up 
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to four ( possibly all) BHs should be eje cted from the clus- 
ter core (Sigurdsso n fc Hernq uist 1993), or similarly that 
nearly all BHs should be ejected, with the possibility of re- 
maining BHs capturing normal stars t o form low-mass X-ra y 
binaries in low-dens ity environments (|Kulkarni et al.lll993l h 
IStrader et all (|2012T ) claim their observations to be in con- 
trast to such prior theoretical predictions. However it was 
not possible to test this in direct iV-body models compa- 
rable to clusters like M22 until very recently as such mod- 
els are computationally expensive. We note that multiple 
BHs have already been found to rema in in iV-body mod- 
els w ith smaller numbers of stars (e.g. iMackev et al.l [20071 . 
120081 ) as well as in Monte Carlos simulations of globula r 
cluster evolution (e.g. iDowning et al.ll2010l : lDowninell2012l h 
wh ile BHs are eje c ted fr om the cluster in iV-body models 
bv lBaneriee et al.1 (|2010h where high numbers of BHs were 
added initially. 

Based on the recent findings by iRepetto et al.l (|2012T l 
suggesting that BHs receive similar velocity kicks as neu- 
tron stars (NSs) upon formation, and the observations 
that not all NS s are expected to receive a high kick 
jPfahl et alj|2002h . we present a direct iV-body model for an 
intermediate-mass globular cluster containing 262 500 stars 
in total and incorporating a retention fraction of w 10% 
for both BHs and N Ss. The cluster is evolved using NB0DY6 
(|Aarsethlll9"99l . 120031 ') in a Milky Way-like gravitational po- 
tential field. 

While the model presented in this paper is by no means 
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an attempt at a direct model of M22, it is representative ow- 
ing to comparable dynamical and absolute ages. We evolve 
the model up to an age of 20 Gyr, and focus our analysis at 
the age of 12 Gyr: the estimated age for many globular clus - 
ters in the Milky Way including M22 i|Salaris fc Weissll2002h . 
At this stage, « 180 000 stars are still retained within the 
cluster (including 16 BHs), implying that this is the largest 
dir ect iV-body model of a globular cluster currently ava ilable 
(cf. iBaumeardt fc Makinolliool : iHurlev fc Sharall2012h . 



2 SIMULATION METHOD & CHOICE OF 
PARAMETERS 

We use the direct TV-body code NBDDY6 ( | Aarsethll 199fil 120031 ; 
iNitadori fc Aarsethl |2012| ) to evolve our mod el and chose 
a set- up similar to the clusters presented in ISippel et al.l 
121)121. We evolve a star cluster with AT, = 250 000 stel- 
lar systems and an initial binary fraction of bf = 0.05, 
implying that the initial numbers of stars is in fact N = 
262 500 including 25 000 stars in 12 500 binary systems. The 
stars are drawn from the initial mass function presented by 
iKroupa et all (|l993h within the limits 0.1 < m < 50 M 
and distributed according to a Plu mmer density profile 
jPlummerl Il91ll ; lAarseth et "all Il974 ). The initial cluster 
mass adds up to Mi — 1.6 x 10 s Mq . The conversion from N- 
body units to physical scale sizes leaves the scale radius as a 
free parameter which is set to give an initial half-mass radius 
of r 50% = 6.2 pc. The cluster is orbiting a Milky Way like 
galactic potential consisting of three c omponents: a point- 
mass bulge, an extended smooth disc l|Mivamoto fc Nagail 
Il975l ) and a logarithmic dark matter halo (Aarscth 2003). 
We chose a gal actocentric distance of d gc = 8.5 kpc. From 
King model fits (Kins 1966) we find an initial tidal radius of 
rt ~ 50 pc for the iV-body model, while M22 has a smaller 
tidal radius of r t = 27 pc |Mackev fc van den BerghlliooBr i 
owing to the smaller galactocentric distance of d gc — 4.9 pc. 

Within the Milky Way, globular clusters across the 
whole range of metallicities exist around the galactocen- 
tric distance used here and we chose a low to intermedi- 
ate metallicity Z= 0.001 (corresponding to [Fe/H]~ —1.3). 
This metallicity is close to that of M22 [Fe/H]= -1.7 
and stars within this metallicity r ange evolve in an almost 
identical fashion (as illustrated in lHurlev et alJfeOQCl . 120041 ; 
ISippel et alj2012T ). Stars are e volved according to stellar and 
binary evolution algorithms (|Hurlev et al. 2000, 2002) im- 
plying that stellar collisions may occur and stellar remnants 
such as white dwarfs (WDs), NSs and BHs form. BHs form 
from progenitors with masses greater than 20 Mq, where the 
resulting remnants are within the mass ran ge of 4 — 30 Mq 
jHurlev et all |2000| ; iBelczvnski et afll201Ch . These massive 
stars evolve rapidly with most BHs forming during the first 
« lOMyr of cluster evolution. 

2.1 Velocity kicks 

It is well established that some NSs will receive velocity 
kicks upon formation re sulting from asym metries in the stars 
surface before collapse (|Pfahl et al.ll2002T ). To mimic a neu- 
tron star retention fracti on similar to indi c ations from ob- 
servations of w 10% fc.f. iPfahl et~ai1l2002l ; |PfahJl2003h we 
adopt a random fiat kick distribution between — 100 km/s. 




Figure 1. Location of BHs (red crosses) at 12 Gyr of cluster age, 
projected in front of the cluster to increase visibility. 

Recent indications suggest that BHs should receive simi- 
lar kicks to NSs (jRepetto et al. 1 120121 ). as such we apply 
an identical kick distributio n to BHs. F or our model, the 
escape velocity is v csc = y /, 2GM/r 50 % = 14 km/s at the 
very beginning of cluster evolution (M is the cluster mass, 
r 50% the half-mass radius and G the gravitational constant). 
At an age of 200 Myr, by which time all BHs and NSs will 
have safely formed, the escape velocity is 10.9 km/s owing 
to r 50% = 8.5 pc and M = 11.8 x 10 4 Mq. 

We note that in general, for a supernova occurring in 
a tight binary system, the kick velocity can be significantly 
above the cluster escape velocity and still result in one or 
both binary components being retained, with the binary 
loosened or broken up. Seven BHs form while part of a bi- 
nary system, all of them being disrupted either immediately 
or within the next « 10 Myr. Ultimately, this means that 
out of « 450 stellar remnant BHs formed through stellar 
evolution, 48 are retained beyond the first 200 Myr, i.e. are 
not immediately ejected from the cluster via the velocity 
kick process. Within less than one Gyr, they have settled 
into the cluster core owing to mass-segregation and have a 
steadily decreasing half-mass radius of ~ 1 pc. The spatial 
distribution of BHs within the model cluster at 12 Gyr is 
illustrated in Fig. [T] (see also Fig. [2}. 

2.2 Size and time scales: measurements 

The sizes of iV-body star cluster models are usually mea- 
sured by either the half-mass radius or the core radius. Both 
quantities are measured in three dimensions, in contrast to 
observational size scales. The half-mass radius r 50 % is sim- 
ply the 50% Lagrangian radius, however the iV-body core 
radius is a density weighted average distance from the clus- 
ter centre ijvon Hoernerlll96ol . 1 19631 ; ICasertano fc HutJll985l ; 
lAarseth! 120031 ') and hence not comparable to the observa- 
tional core radius resulting from e.g. King model fits to the 
surface brightness profile. The procedure to calculate both 
observable half-light as well as core radii is similar to that 
described in lSippel et al.l (|2012h . In particular, we cross con- 
volve the NB0DY6 output of mass, luminosity and radius for 
each single st ar with stellar atmosphere model calculations 
l|Kuruczll979h to obtain V-band magnitudes. Surface bright- 
ness profiles are produced by taking the average over pro- 
j ection along the x—, y — and z— axes; and we use gridfit 
l|McLaughlin et alj [2008'l to fit a King model to the cluster. 
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Figure 2. Size evolution of the cluster over time, from top to 
bottom these are half-mass r 5 o% (dashed), TV-body core radius 
r c (dotted), further Lagrangian radii r x %, r 5% and r 1% as 
well as the half- mass radius for the BH population (red), always 
smaller than even the 0.1% Lagrangian radius. 

From projection effects alone, the half-mass radius is ~ 1.3 
times larger than the half-light radius. 



3 EVOLUTION 

The first few Gyr of cluster evolution are dominated by 
heavy mass loss arising from stellar evolution as well as dy- 
namical relaxation that results in an initial expansion of the 
cluster to a half-mass radius of up to r 50 % ~ 11 pc (Fig. 
[2J. The cluster then slowly starts to contract to r 50 % = 
10.6 pc and an TV-body core radius r tc = 3.1 pc, both three- 
dimensional quantities, at 12 Gyr. In terms of physical size, 
the model has a half light radius of Th = 6.8 pc at 12 Gyr 
and we find an observational core radius of r oc = 3.6 pc. The 
half- mass relaxation time at 12 Gyr is « 2.1 Gyr, where we 
utilize the h alf-light rad i us to to c alculate this n umber as 
suggested by iDiorgovskil (|l993l ) and iHarrisl l|l996l ). 

Owing to the influence of the external galactic tidal field 
as well as dynamical interactions, stars are continuously lost 
from the cluster, with ~ 180 000 stars remaining at 12 Gyr. 
This corresponds to w 68% with respect to the initial num- 
ber of stars and a total mass of M = 6.7 x 10 4 M Q (« 42% 
of the initial mass). 

Globular clusters are efficient in equalizing the distribu- 
tion of energy, causing the most massive stars or remnants 
to sink to the centre, resulting in mass segregation. Even 
though our model is segregated after 12 Gyr, the average 
mass within the innermost pc is only 0.7 Mq (compared to 
an overall average of 0.38 Mq). Simultaneously, heating of 
the core from either BHs or NSs can be expected to produce 
signa tures such as an expanded core (jMackev et al.l 120071 . 
l200Sf ) . where both effects can influence the surface bright- 
ness profile and are measurable by means of the core radius 
l|Kindll966l ). M22 has a large core radius compared to other 
globular clusters in the Milky Way, indicating that a heating 
source may be present. While our model does not enter the 
phase of core collapse during its evolution, the subsystem 
of stars contained within the innermost 2pc (dynamically 
dominated by BHs) is entering a phase of core-collapse at 
~ 1 Gyr of cluster age, indicated by the drop and fluctua- 
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Figure 3. Contribution of stars to the central parsec of the clus- 
ter, i.e. d < 1 pc, where the distance d is again measured in three 
dimensions. It can be seen that at 12 Gyr, MS stars as well as 
WDs and BHs are contributing roughly 1/3 of the dynamical 
mass each, while the contribution of other stars (giants or neu- 
tron stars) is minimal (i.e. < 10 Mq for both together). The ratio 
is different when looking at the number of stars in each category: 
at 12 Gyr, MS stars count for roughly 2/3 of the stars within the 
innermost pc, while WDs make up fa 1/3 of the stars. The dip in 
BH mass within the central region at 5.5 Gyr is caused by an 
energetic interaction within a close subsystem of three BHs, two 
of them previously in a BH-BH binary. One BH is ejected out 
of the cluster, while the BH-BH binary receives a recoil velocity, 
ejecting it momentarily beyond the central parsec. The loss of ki- 
netic energy causes an expansion of the cluster core, however the 
BH-BH binary quickly sinks back to the centre. 

tion in the core radius (Fig. [2]). This results from the fact 
that the BH population has a much smaller half-mass re- 
laxation time than the cluster as a whole, resulting the the 
BHs segregating to the cluster centre. The half-mass radius 
of the entire BH population is ~ 0.2 pc and smaller than 
the 0.1% Lagrangian radius (see Fig. [2}. At 12 Gyr, the to- 
tal mass contained in BHs is ~ 190 Mq corresponding to 
only « 0.3% of the overall cluster mass but ~ 32% of the 
mass within the central parsec. 

This already indicates that BHs are an important com- 
ponent to the contribution of mass in the cluster centre. In- 
deed we find that at 12 Gyr, the innermost parsec (measured 
in three dimensions) is composed equally of main sequence 
stars, white dwarfs, and BHs (Fig. [3j) . For the metallicity of 
both M22 and our model the main sequence turnoff mass is 
0.83 Mq at 12 Gyr, implying that most white dwarfs (with 
a peak in their mass distribution at 0.6 Mq) will be less 
massive than the main sequence turnoff, while all neutron 
stars and black holes will be more massive. In total, 245 Mq 
is contained in the innermost pc in 352 stars (compared to 
2 322 stars or 1 263 Mq in projection), and this central den- 
sity increases to 528 stars or 423 Mq at an age of 20 Gyr (or 
2 582 stars and 1616Mq in projection). 

3.1 Binary systems with black holes 

Binary formation and disruption, as well as hardening of bi- 
nary systems, are common processes in globular clusters. 
However it is easier to form new binary systems via ex- 
change interactions than creating them in two-body encoun- 



4 Anna C. Sippel, Jarrod R. Hurley 



ters (|Heggiell 19751 ; iHeggie fc Hudl2003l ). A black hole can be 
part of a binary system where the other component can be 
at any stellar evolution stage. In our simulation, the first dy- 
namical binary system including a BH forms at ~ 380 Myr 
(a 25.3 M @ BH and a 0.2 Mq MS star). At I2Gyr, the ini- 
tial binary fraction of bf = 0.05 is slightly decreased to 
bf = 0.045 (with 7850 binary systems present). In M22, 
one or potentially both observed BHs are in a binary with 
a low-mass main se quence star, or possibly a white dwarf 
jStrader et al.ll2012h . In our model, we find that at around 
12 Gyr, two BHs are in a long-period binary system with a 
main sequence star each, where neither of those systems are 
undergoing mass transfer. Specifically, system a) consists of 
a BH with mass 6.2 Mq and a main sequence star with mass 
0.8 Mq while system b) consists of a BH with mass 6.4Mg 
and a main sequence star of 0.4 Mq . Just before 12 Gyr a 
three-body encounter occurs involving system b) and a main 
sequence star with mass 0.7 Mq, replacing the lighter main 
sequence star and forming again a binary system consisting 
of a black hole and main sequence star. 

A BH-BH binary system is also present at around 
12 Gyr, where the components of this binary have masses 
of 22 Mq each. Just after 12 Gyr, during a three-body en- 
counter with another BH with 13 Mq , this binary gets dis- 
rupted, whilst the individual components of the three-body 
system stay retained within the cluster. A BH-BH binary 
in a new combination forms quickly and the total number 
of 16 BHs is stable around the 12 Gyr timeframe, as illus- 
trated in Fig. [4] Dynamical interactions like these are what 
can ultimately cause the depletion of all but one BH or a 
tight BH-BH binary (which might ultimately merge, lAarsethl 
|2012| ). from the cluster. However it is not certain over which 
timescale the evaporation of black holes from globular clus- 
ters should take place. This depletion timescale is dependent 
on random encounters, necessary to eject BHs, while the 
frequency of such encount ers depends on clus ter parameters 
such as the core density (|Poolev et alj|2003h . We conclude 
that 12 Gyr is not enough time to deplete all or almost all 
BHs. In fact we have evolved the model up to 20 Gyr, and 
even at this late stage, 10 BHs remain in the cluster. Five 
of those are in a binary system with either a main sequence 
star, white dwarf or another black hole. We also find a BH- 
NS system with a total mass of ~ 18 Mq at ~ 18 Gyr in a 
stable configuration for « 0.6 Gyr, at which stage the NS 
gets replac ed by a BH w i th m ass 7.1 M© in a three-body 
encounter. IClausen et al.l (|2012l ) have recently presented a 
study regarding BH-NS mergers in GCs and we conclude 
that dynamical interactions don't efficiently produce BH-NS 
binaries in GCs. 

3.2 Comparison 

While a direct model approach for M22 is still out of reach 
from a computational point of view, the TV-body model pre- 
sented in this work is comparable in many aspects. While 
our model has a half-mass relaxation time of t r h = 2.1 Gyr, 
the correspinding value is i r h = 1.7 Gyr for M22 (|Diorgovskil 
ll993l ; lHarrislll996l ). M22 has a half-light r h = 3.1 and core 
radius r oc = 1.2 pc (|Harris!l 19961 ). implying that the TV-body 
model is less dense than M22. We argue that the higher cen- 
tral density of M22 does not severly change our predictions 
for stellar remnant BHs in globular clusters with large core 
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Figure 4. Evolution of the mass of the entire black hole pop- 
ulation a) and the number of BHs within the cluster b). In the 
upper panel, the solid black line denotes the maximum mass of 
the BHs still retained in the cluster, while the blue-dashed line 
denotes the mean mass of the BH population at a given time. The 
number of BHs is decreasing b) from 48 at 200 Myr to 16 BHs at 
12 Gyr of cluster age. A yellow background denotes the presence 
of one BH-BH binary, while orange means two such systems are 
present, and white none. Note that the disruption of both BH- 
BH binary systems at « 6 and m 11 Gyr is accompanied by the 
ejection of the most massive black hole contained in the cluster. 
The disruption of one BH-BH binary is also accompanied by the 
ejection of one BH at £d 2 and ss 7.5 Gyr. 



radii since examples of denser model clusters with BHs exist 
in the literature. 

Most pro m inent is the TV-body model by 
iHurlev fc Sharal (|2012T ) consisting of 200 000 stars ini- 
tially, evolved at a galactocentric distance of 3.9 kpc and 
with an initial half- mass radius of 4.7 pc. The strong 
gravitational field of the host galaxy at this distance has an 
impact on the cluster by stripping off stars in the outskirts, 
resulting in only ~ 35 000 stars remaining at 12 Gyr. How- 
ever, the half-mass relaxation time of the IHurlev fc Sharal 
model is much below that of M22 resulting in the 
model reaching the end of its initial phase of core-collapse 
at 11.5 Gyr, greatly exceeding the dynamical age of M22. 
At this stage, the core radius is r oc = 0.84 pc with four 
BHs still present (two of those in a BH-BH binary). This 
is the case even though the central density is significantly 
enhanced at this time (just over 3 000 stars in the innermost 
parsec, or « 5 600 in projection). Even though the black 
holes ultimately get ejected during the evolution subsequent 
to core-collapse, we can safely conclude that BHs in the 
core of M22 can be expected , as th e density of M22 lies in 
between the IHurlev fc Sharal i|2012l ) model and the TV-body 
model presented in this work - however much closer to the 
new model here. 

Other examples of TV-body cluster models with BHs 
currently exist in the li terature (jMackev et al.l |2007| . 120081 ; 
IHurlev fc Mackevl l20ld just to name a few), as well as 
Monte Carlo models of more massive globular clusters 
at various densities (|Downing et al ] |2010| ; iDownind |2012l ; 
iMorscher et al.ll2012h . In addition we note that models in- 
volving mergers of stellar remnant black holes as well as well 
as intermediate mass BHs have been carried out in the past 
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(e.g. iPortegies Zwart &: McMillan 



l2004l ; IO'Learv et alj|2006l : lAarsethj f 



.20021 ; iBaumgardt et al.l 
2012T) . 



4 CONCLUSIONS 



We conclude that all three findings of IStrader et alj (|2012l ): 
a) the observation of two stellar remnant black holes in M22 
which are b) potentially in binary systems with main se- 
quence stars and that c) 5 — 100 stellar remnant black holes 
might be present in a cluster such as M22 - are in excellent 
agreement with our direct iV-body model. A key point in 
this analysis is that M22 has a large core radius compared to 
the average for the globular cluster population of the Milky 
Way, and is not dynamically old. We conclude that multi- 
ple stellar remnant black holes in the core of such clusters 
can exist even beyond the current time. Furthermore, in the 
JV-body model presented here we find that 1/3 of the BHs 
retained in the cluster upon formation remain in the core 
beyond the Hubble time. 

With further advances in computing hardware we ex- 
pect these results to be confirmed by even larger models 
with an expanded set of initial parameters in the near fu- 
ture. 
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